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Abstract 

Experimental X-ray crystallography, NMR (Nuclear Magnetic Resonance) 
spectroscopy, dual polarization interferometry, etc are indeed very pow- 
erful tools to determine the 3 -Dimensional structure of a protein (includ- 
ing the membrane protein); theoretical mathematical and physical computa- 
tional approaches can also allow us to obtain a description of the protein 
3D structure at a submicroscopic level for some unstable, noncrystalline 
and insoluble proteins. X-ray crystallography finds the X-ray final struc- 
ture of a protein, which usually need refinements using theoretical proto- 
cols in order to produce a better structure. This means theoretical meth- 
ods are also important in determinations of protein structures. Optimiza- 
tion is always needed in the computer-aided drug design, structure-based 
drug design, molecular dynamics, and quantum and molecular mechanics. 
This paper introduces some optimization algorithms used in these research 
fields and presents a new theoretical computational method - an improved 
LBFGS Quasi-Newtonian mathematical optimization method - to produce 
3D structures of Prion AGAAAAGA amyloid fibrils (which are unstable, 
noncrystalline and insoluble), from the potential energy minimization point 
of view. Because the NMR or X-ray structure of the hydrophobic region 
AGAAAAGA of prion proteins has not yet been determined, the model con- 
structed by this paper can be used as a reference for experimental studies on 
this region, and may be useful in furthering the goals of medicinal chemistry 
in this field. 

Keywords: Protein 3D Structure; Computational Approaches; Optimization Method; 
Molecular Modelling; Prion AGAAAAGA Amyloid Fibrils. 
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Highlights: ► Determinations of 3D structures of unstable, noncrystalline and in- 
soluble proteins. ► Theoretical computer, mathematical/physical computational 
approaches and concepts. ► Time-consuming, costly X-ray crystallography and 
NMR spectroscopy powerful tools. ► Mathematical local search optimization 
methods to solve three-body physical problems. ► Using computational chem- 
istry Amber/Gromacs packages to refine molecular models. 

1 Introduction 

Neurodegenerative diseases including Parkinson's, Alzheimer's, Huntington's, and 
Prion's were found they all featured amyloid fibrils [HJ |2j |3] ffl El O . Amyloid is 
characterized by a cross-/? sheet quaternary structure and recent X-ray diffraction 
studies of microcrystals revealed atomistic details of core region of amyloid [0[8). 
All the quaternary structures of amyloid cross-/? spines can be reduced to one of 
the 8 classes of steric zippers of [8 |, with strong van der Waals (vdW) interactions 
between /3-sheets and hydrogen bonds (HBs) to maintain the /3-strands. A new 
era in the structural analysis of amyloids started from the 'steric zipper'- /3-sheets 
(7]|. As the two /3-sheets zip up, Hydrophobic Packings (HPs & vdWs) have been 
formed. The extension of the 'steric zipper' above and below (i.e. the /3-strands) 
is maintained by Hydrogen Bonds (HBs) (but usually there is no HB between the 
two /3-sheets). This is the common structure associated with some 20 neurode- 
generative amyloid diseases. 

We first do some mathematical analysis for the common structure. Let r be 
the distance between two atoms, the vdW contacts of the two atoms are described 
by the Lennard- Jones (LJ) potential energy: 

VUr)=4e((-) 12 -(-) 6 ), (1.1) 
\ r r J 

where e is the depth of the potential well and o is the atom diameter; these pa- 
rameters can be fitted to reproduce experimental data or deduced from results of 
accurate quantum chemistry calculations. The (^) 12 item describes repulsion and 
the —(-) 6 item describes attraction (Figure 1). If we introduce the coordinates of 
the atoms whose number is denoted by iV and let e = a = 1 be the reduced units, 
then, for iV atoms, Eq. 1.1 is 

- - / 1 1 \ 

/(*) = *£ £ (□r-isj. ( L2 ) 

i=i j=i,j<i \ H y/ 

where = (x 3 i^ 2 -x 3j ^2) 2 +{x 3 i-i-x 3j -i) 2 +(x 3 i-x 3j ) 2 , and (x 3i _ 2 , x 3i _i, x 3i ) 
is the coordinates of atom i, N > 2. The minimization of LJ potential energy f(x) 
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Figure 1: The LJ potential energy (Eqs. 1.1 and 1.4) (can be seen in Figure 1 of 



on W n (where n = 3N) is an optimization problem 



min fix) 



;i.3) 



which is a well-known and challenging test problem for global optimization (see 
|http://www-wales.ch.cam.ac.uk/CCD.htmI1 and its recent references such as [10] QT] 
d2l[T3l[T4l). It is very hard for global optimization to directly solve Eq. 1.3 even 
without a large number of atoms. Similarly as Eq. 1.1 - i.e. 
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the potential energy for the HBs between /3-strands is 



V HB (r) 



C_ 

Zl2 



D 

^10' 



1.5) 



where A, B, C, D are given constants and usually most of the HBs are still kept 
during the phase of molecular modeling. Thus, the amyloid fibril molecular mod- 
eling problem can be reduced to solve the optimization problem Eq. 1 .3 though it 
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is not easy to accurately solve Eq. 1 .3 for a large molecule. 

Alternatively, we have found another way to solve Eq. 1.3 |[T5l . Seeing Figure 
1 , we may know that the optimization problem Eq. 1 .3 reaches its optimal value at 
the bottom of the LJ potential well, where the distance between two atoms equals 
to the sum of vdW radii of the two atoms. Hence, the amyloid fibril molecu- 
lar modeling problem can be looked as a molecular distance geometry problem 
(MDGP) |[T6ll . As an example to explain MDGP, the problem of locating sensors 
in telecommunication networks is a DGP. In such a case, the positions of some 
sensors are known (which are called anchors) and some of the distances between 
sensors (which can or cannot be anchors) are known. The DGP is to locate the 
positions of all the sensors. The MDGP looks sensors as atoms and their telecom- 
munication network as a molecule. In mathematics, the following Eqs. 1.6~1 .8 
can express the MDGP for Eq. 1.3. The 3D structure of a molecule with N atoms 
can be described by specifying the 3D coordinate positions x\, x 2 , ■ ■ • , xjy G M 3 
of all its atoms. Given bond lengths between a subset S of the atom pairs, the 
determination of the molecular 3D structure is 

(Vq) to find xi,X2,---,xn such that \\xi — Xj\\ = dij, (i, j) G S, 

(1.6) 

where 1 1 • 1 1 denotes a norm in a real vector space and in this paper it is calculated 
as the Euclidean distance 2-norm. Eq. 1.6 can be reformulated as a mathematical 
global optimization problem 

(V) mmP(X) = ^i(lki-^H 2 -4) 2 (1-7) 

in the terms of finding the global minimum of the function P(X), where Wij, G 
S are positive weights, X = (xi,x 2 , . . . ,xn) t G W 1 [fTTl and usually S has 
fewer elements than N 2 /2 due to the error in the theoretical or experimental data 
|[T8l [TBI . Even there may not exist any solution xi, x 2 , ■ ■ ■ , x^ to satisfy the dis- 
tance constraints in Eq. 1.6, for example when data for atoms i,j, k G S violate 
the triangle inequality; in this case, we may add a perturbation item — e T X to 
P(X): 

(V e ) mmP £ (X)= w ij {\\x i -x j \\ 2 -d%f-e T X : (1.8) 
(i,j)es 

where e > 0. In some cases, instead exact values d^, G S can be found, 
we can only specify lower and upper bounds on the distances: Uj < \ \xi — Xj\ | < 
Uij, G S; in such cases we may penalize all the unsatisfied constraints into 
the objective function of (V £ ) by adding J2a,j)eS ( max — — x jll 2 50}) 2 + 
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(max { 1 1 x i — Xj 1 1 2 — uf 3 ; , } ) into P £ (X) OH [16l , where we may let dij be the 
interatomic distance (less than 6 angstroms) for the pair in successive residues of 
a protein and set = (1 — 0.05)^ and Uij = (1 + 0.05)^ IfToTl . In this paper, 
we aim to solve Eq. 1.8 (or Eq. 1.3) for modeling amyloid fibril molecular 3D 
structures. 

Neurodegeneration is the progressive loss of structure or function of neurons, 
including death of neurons. A prion is a misshapen protein that acts like an in- 
fectious agent (but not requiring either DNA, RNA, or both) to cause a number 
of fatal diseases. Prion diseases are rich in /3-sheets (compared with the normal 
prion protein PrP c in rich of a-helices) and are so-called "protein structural con- 
formational" diseases. The normal hydrophobic region 113-120 AGAAAAGA 
peptide of prion proteins is an inhibitor/blocker of prion diseases. PrP lack- 
ing this palindrome could not convert to prion diseases. Brown et al. pointed 
out that the AGAAAAGA peptide was found to be necessary (though not suf- 
ficient) for blocking the toxicity and amyloidogenicity of PrP 106-126, and the 
peptide AGAA does not form fibrils [[T9l . The minimum sequence necessary for 
fibril formation should be AGAAA, AGAAAA, AGAAAAG, AGAAAAGA and 
GAAAAGA, but the molecular structures of these fibrils have not known yet. This 
paper addresses an important problem on modeling the 3D molecular structures 
of prion AGAAAAGA amyloid fibrils of neurodegenerative diseases. The rest of 
this paper is arranged as follows. In the next section, i.e. Section 2, an improved 
LBFGS Quasi-Newtonian method is presented for solving Eq. 1.8. Section 3 im- 
plements this Quasi-Newtonian method by constructing an 3D molecular structure 
of prion AGAAAAGA amyloid fibrils of neurodegenerative prion diseases. Nu- 
merical results of computations show that the method designed in Section 2 is 
very effective and successful. This concluding remark will be made in the last 
section, i.e. Section 4. 

2 Methods 

In a (macro)molecular system, if it is very far from equilibrium, then the forces 
may be excessively large, a robust energy minimization (EM) is required; another 
reason to perform an EM is the removal of all kinetic energy from the system: 
EM reduces the thermal noise in the structures and potential energies [|20l . EM, 
with the images at the endpoints fixed in space, of the total system energy pro- 
vides a minimum energy path. EM can be done using steepest descent (SD), con- 
jugate gradient (CG), and Limited-memory Broyden-Fletcher-Goldfarb-Shanno 
(LBFGS) methods. 
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Three kinds of possible EM methods are: (1) derivative-free methods - that 
require only function evaluations, e.g. the simplex method and its variants; (2) 
derivative information methods - the partial derivatives of the potential energy 
with respect to all coordinates are known and the forces are minimized, e.g. 
SD, CG methods; and (3) second derivative information methods, e.g. LBFGS 
method. "SD is based on the observation that if the real- valued function f(x) is 
defined and differentiable in a neighborhood of a point x then f(x) decreases 
fastest if one goes from x in the direction of the negative gradient of f(x) at xq 
and SD local search method converges fast lETI . SD is robust and easy to im- 
plement but it is not most efficient especially when closer to minimum; at this 
moment, we may use the efficient CG. CG is slower than SD in the early stages 
but more efficient when closer to minimum. CG algorithm adds an orthogonal 
vector to the current direction of the search, and then moves them in another di- 
rection nearly perpendicular to this vector. The hybrid of SD-CG will make SD 
or CG more efficient than SD or CG alone. However, CG cannot be used to find 
the EM path, for example, when "forces are truncated according to the tangent 
direction, making it impossible to define a Lagrangian" [B21I231 . In this case, the 
powerful and faster quasi-Newtonian method (e.g. the LBFGS quasi-Newtonian 
minimizer) can be used E2 |24l 1251 |26l |27l |28]|. We briefly introduce the LBFGS 
quasi-Newtonian method as follows. 

Newton's method in optimization explicitly calculates the Hessian matrix of 
the second-order derivatives of the objective function and the reverse of the Hes- 
sian matrix [|29l . The convergence of this method is quadratic, so it is faster than 
SD or CG. In high dimensions, finding the inverse of the Hessian is very expen- 
sive. In some cases, the Hessian is a non-invertible matrix, and furthermore in 
some cases, the Hessian is symmetric indefinite. Quasi-Newton methods thus ap- 
pear to overcome all these shortcomings. 

Quasi-Newton methods (a special case of variable metric methods) are to ap- 
proximate the Hessian. Currently, the most common quasi-Newton algorithms 
are the SRI formula, the BHHH method, the widespread BFGS method and its 
limited/low-memory extension LBFGS, DFP, MS, and Broyden's methods |[30l 
[HI [32l[3l. In Amber [J23]| and Gromacs [HI, LBFGS is used, and the hybrid of 
LBFGS with CG - a Truncated Newton linear CG method with optional LBFGS 
Preconditioning [25 J - is used in Amber [1231 . 

For BFGS method, whether it converges at all on nonconvex problems is still 
an open problem. In fact, Powell (1984) gave a counter-example that shows that 
BFGS with an inexact line-search search may fail to converge 041 l35l l36l . Li 
and Fukushima (2001) proposed a modified BFGS method for nonconvex objec- 
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tive function [T3VTI . Basing on 11371 1241 [38l 1251 [39ll , in this paper we present an 
improved LBFGS method described as follows fl40l - which presents the non- 
monotone line search technique [HT1I421 for the Wolfe-type search. The improved 
LBFGS method presented in this paper is much better than the standard BFGS 
method in view of the CPU time (see Figure [2]) tested through more than 30 non- 
linear programming problems (where each selected problem is regular, that is, its 
first and second derivatives exist and are continuous everywhere, and each prob- 
lem is with different dimensions, i.e., 100, 500, 1000 and 10000 dimensions) and 
its mathematical theory to support this algorithm can be seen from the Supple- 
mentary Materials 11401 listed at the end of this paper. This paper implements the 
Wolfe-type search by the approximation technique of piecewise linear/quadratic 
function 11431 . 

Algorithm 1 An Improved LBFGS Method for minimizing nonconvex func- 
tion gol 

Step 0: Choose an initial point xq G IR n , an initial positive definite matrix H , 
and choose constants <J\,<j% such that < <j\ < a 2 < 1, and choose an positive 
integer ra\. Let k = 0. 

Step 1: If ||<7fc|| = 0, then output x k and stop; otherwise, go to Step 2. 
Step 2: Solve the following linear equation to get d k : 

dk = —Hk9k, 

Step 3: Find a step-size \k > satisfying the Wolfe-type line search conditions: 

f(%k + A fc 4) < f(xk) + criAfc#fc 4, (2.1) 

g(x k + \ k d k ) T d k > a 2 g T k d k . (2.2) 
Step 4: Let the next iterate by x k+ \ = x k + X k d k . Calculate g k+ i and \\g k+ i\\. 
Step 5: Lets k = x k+1 -x k = X k d k ,y k = g k+ \-g k , 7fc = ||S7c||> theny* k = y k +^ k s k . 
Step 6: Let m = mm{k + 1, mi}. Update H k following the formula 

H k+ i = (V fc * T ---V fc % l+ Jif fc _m+i(V fc t_^ +1 ---V fc *) 

+ W fc-m+l(^A*-l ' ' ' ^ k ^m+2) S k-rh+lS k _ ih+1 {y^ m+2 ■ • • V^*_l) 

+ • • • + u* k s k sl, 

(2.3) 

where u* k = V k * = I - u* k y* k s T k (when k = 0, H = V* = I - 

44, Hl = (I- k^) T H (I - 44) + ^s k s T k ). 
Step 7: Let k = k + 1 and go to Step 1. 
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In Algorithm 1, g k denotes the gradient of f(x) at x k , and the convergence and 
the R-linear convergent rate of Algorithm 1 are guaranteed by the following as- 
sumptions: (Al) / is twice continuously differentiable, (A2) / has Lipschitz con- 
tinuous gradients and Hessians, (A3) the Hessian at the stationary point is always 
positively definite, and (A4) the level set of / is bounded. The detailed proof of 
the convergence and i?-linear convergent rate of Algorithm 1 can be found in the 
Supplementary Materials supplied at the end of this paper. 




Value ol t 



Figure 2: Performance based on CPU time, where the definitions of P and r can 
be found in 144]. 



In BOl . the following nonmonotone line search technique for stepsize \ k is 
used: 

f(x k + X k d k ) < max f{x k -j) + oi\ k gld k , 

0<j<Mo 

g(xk + \d k ) T d k > max{cr 2 , 1 - (A*| \d k \ \) p }gld k , 

where p E (— oo, 1), and M is a nonnegative integer. This is a difference between 
the algorithm of ll40l and this paper. All in all, it is well known that quasi-Newton 
method is an efficient solution method for unconstrained and continuously dif- 
ferentiable minimization problem B31l46ll47l . However, it needs computing and 
storage of the updated matrix which is an approximation to the Hessian matrix at 
each iteration of the method. Hence, its efficiency may decrease when it is applied 
to large scale optimization problem. To overcome the drawback, limited memory 
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quasi-Newton method is proposed [|48Tl . The main ideal of this method is nearly 
identical to that of the standard BFGS method, and the only difference is that the 
inverse Hessian approximation is not formed explicitly, but defined by a small 
number, say in, of BFGS updates. This technique received much attention in re- 
cent years and numerical experiments show that it is very competitive ll49l l24ll . 
and its global convergence and R-linear convergence rate with Wolfe line search 
are established for the uniformly convex case 11261 1241 . Since the limited mem- 
ory BFGS method may suffer from ill-conditions for small value of m, Al-Baali 
(1999) ll50l made some modifications to the method and establish its global con- 
vergence based on the same assumptions, and Byrd et al. (1994) [51] derives new 
representation of limited memory quasi-Newton matrices for the benefit of com- 
puting the updated matrix. Recently, a non-monotone line search is introduced, 
see e.g., ||4T1I4"21 . Then it is showed to be more competitive and practical for solv- 
ing nonlinear optimization problems, and [1521 established the global convergence 
of this line search applied to limited memory BFGS method based on the uni- 
formly convex assumption. Motivated by the above observation, it turn out that in 
two respects the limited memory BFGS method is much less effective. First, we 
note that the convergence analysis of these method are focused on the uniformly 
convex assumption and little is known for nonconvex case. Second, numerical 
experiments have suggested the main weakness of limited memory method is that 
it may converge very slowly in terms of number of iterations for ill-conditioned 
problems. The purpose of the above Algorithm 1 is to reduce these defects and 
Figure 2 shows the effectiveness of the proposed algorithm. We will apply it into 
the molecular modeling of prion AGAAAAGA amyloid fibrils in the next section. 



3 Results and Discussion 

From their research of prion, scientists found that the cross-/? structure of peptides 
is with the nature of self-aggregation, the self-aggregating to form fibers. This 
provides us a new research idea for nanomaterials. HBs can be formed between 
peptide /3-strands, and one peptide monomer connects together with another in ac- 
cordance with the specific structure to form fibers. Many laboratories in the world 
are synthesizing peptides that can self-aggregate to form fibers, and want to be 
able to control the growth of the fiber to find out new functional materials [|53ll54"l . 
The studies of this paper not only benefit nanometerials research, but also benefit 
the research on neurodegenerative amyloid fibril diseases. Prion AGAAAAGA 
peptide has been reported to own an amyloid fibril forming property (initially de- 
scribed in 1992 by Gasset et al. of Prusiner's Group) [fl9l 1551 1551 1571 [581 1591 [601 
EB M M M M M M M MMUIlMn^Mn^MniinM, but there has 
not been experimental structural bioinformatics for this segment yet due to the 
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unstable, noncrystalline and insoluble nature of this region. Furthermore, Zhang 
(201 1) did accurate calculations to confirm the amyloid fibril property at this re- 
gion (Figure 3) [1791 . 



Amyloid fibril identification for prions 



i i i i i i i i i i i r 



human prion 

mouse prion 

rabbit prion 




20 40 60 80 100 120 140 160 180 200 220 240 260 

Residue Number 



Figure 3: Prion AGAAAAGA (113-120) segment is clearly and surely identified 
as the amyloid fibril formation region, because its energy is less than the amyloid 
fibril formation threshold energy -26 kcal/mol [80]. 



3.1 Material for the Molecular Modeling 

This paper uses a suitable pdb file template 20MP.pdb (the LYQLEN peptide 
derived from human insulin residues 13-18 [8|) from the Protein Data Bank to 
build an 8-chain AGAAAAGA prion amyloid fibril molecular model to illuminate 
Algorithm 1 works very well. To choose 20MP.pdb (Figure 4) as the modeling 
template is due to it can pass all the long procedures of SDCG-SA (equilibrations 
& productions)-SDCG of 11791 . By observations of Figure 4 and the 2nd column 
of coordinates of 20MP.pdb, we know that E(F) chains can be calculated on the 
XZ-plane from A(B) chains by Eq. 3.1 and other chains can be got by a parallel 
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beta -sheet 1 



beta -sheet 2 



Figure 4: Protein fibril structure of human insulin LYQLEN (13-18) (PDB id: 
20MP). The dashed lines denote the HBs between the pairs of f3 -strands. A, 
B, C, D, E, F, G, H denote the chains of the fibril. The pair of -sheets 1 & 2 
forms a completely dry interface by vdWs, and between many pairs of (3-sheets 
wet interfaces are formed with water molecules. 

up (or down) along the X-axis by Eqs. 3.2~3.3: 



3.2 New Molecular Modeling Homology Model 

Basing on the template 20MP.pdb from the Protein Data Bank (Figure 4), the 
AGAAAAGA palindrome amyloid fibril model of prions (denoted as Model 1) 
will be constructed. Chains AB of Model 1 will be got from AB Chains of 
20MP.pdb using the mutate module of the free package Swiss-Pdb Viewer (SPDBV 
Version 4.01) (http://spdbv.vital-it.ch). It is pleasant to see that some HBs are still 
kept after the mutations; thus we just need to consider the vdWs only. Making 
mutations for EF Chains of 20MP.pdb, we can get the EF Chains of Model 1. 
Then we add GLY and ALA residues by XLEaP module of Amber 1 1 . However, 
the vdWs between Chain A and Chain E, between B Chain and F Chain are too 



E{F) = A(B) + (-1.885, 0, 17.243) 



(3.1) 



C(D) = A(B) + (9.666,0,0), 
G(H) = E(F) + (9.666,0,0). 



(3.2) 
(3.3) 
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far at this moment (Figure 5, where the twice of the vdW radius of CB atom is 3.4 
angstroms). 



In 11791 the commercial package Insight!! (http://accelrys.com/) is used to build 
models. Instead of InsightH, because this package is not available by the authors, 
this paper uses Algorithm 1 to build and optimize Model 1. In "Zipper 1", fix- 
ing the coordinates of A.ALA3.CB, A.ALA1.CB and letting the coordinates of 
E.ALA6.CB, E ALA4.CB be variables, we may get an optimization problem: 

min 4 | 
+ 4{ 
+4 i 





-1.071) 2 - 


h(xi2-2.986) 2 - 
1 




-1.888) 2 ] H 


l(xu 


-1.071) 2 - 


h(a;i2 


-2.986) 2 +(xi 3 -1.888) 2 ] :i 
1 


[(X21 


-1.071) 2 - 


h(x 22 -2. 986) 2 - 
1 


h(x23 


-1.888) 2 ] 6 


l(X21 


-1.071) 2 - 


<r(.X23 


-2.986) 2 +(a;23-1.888) 2 ] 3 
1 



or 



mini {(xn _ 1.071) 2 + (x 12 - 2.986) 2 + (x 13 - 1.888) 2 - 3.4 2 } 2 
+\ {{x 2l - 1.071) 2 + (x 22 - 2.986) 2 + (x 23 - 1.888) 2 - 3.4 2 } 2 
+i {(x 21 - 1.135)2 + (x 22 + 0.763)2 + (x23 _ 7>209) 2 _ 3>4 2}2 

-0.05 {Xu + Xu + Xl3 + X 2 1 + X 2 2 + X 23 } 



(3.4) 



(3.5) 



with an initial solution (-0.067, 5.274, 7.860; -1.119, 1.311, 13.564). Similarly, 
in "Zipper 2", fixing the coordinates of B.ALA4.CB, B.ALA6.CB and letting the 
coordinates of F.ALA1.CB, F.ALA3.CB be variables, we get another optimization 
problem: 



min 4 



5.446) 2 + (xi2-2.796) 2 +(xi 3 -2.662) 2 ] 6 [(xii-5.446) 2 + (xi2-2.796) 2 +(a'i 3 -2.662) 2 ] 3 

14 J 1 _ 1 

1 [(i2i-5. 446) 2 + (x22-2.796) 2 +(a;23-2.662) 2 ] c, [(x 2 i-5.446) 2 + (x 2 2-2.796) 2 +(a; 2 3-2.662) 2 ] :i 

+4 I 



1 



[(X2i-5.201) 2 + (x22+1.125) 2 +(a;23-7.873) 2 f [(x 2 i-5.201) 2 + (x 2 2+1.125) 2 +(x23-7.873) 2 ] 3 

(3.6) 



or 



min| {(x n - 5.446) 2 + (x 12 - 2.796) 2 + (x l3 - 2.662) 2 - 3.4 2 } 2 
+i {(x 2 i - 5.446)2 + (x 2 2 - 2.796)2 + (x 23 - 2.662) 2 - 3.4 2 } 2 
+i {(£21 - 5.201)2 + (x 22 + 1.125)2 + (x23 - 7.873) 2 - 3.4 2 } 2 
-0.05 {xu + X12 + X13 + X21 + X22 + X23} 



(3.7) 



with an initial solution (4.714, 4.878, 8.881; 4.170, 1.360, 14.292). Next, we 
solve Eqs. 3.5 and 3.7 by Algorithm 1. 

We first solve Eq. 3.5 in the use of Algorithm 1. We set <j\ = 10 -4 , 02 = 
0.1 ~ 0.9, mi = 3 ~ 7, take the initial solution x =(-0.067, 5.274, 7.860; - 
1.119, 1.311, 13.564) and calculate its gradient ^=(-69.7747, 140.135, 365.852, 
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Figure 5: 5a shows the distances of "Zipper 1 "-E.ALA6.CB-A.ALA3.CB- 
E.ALA4.CB-A.ALA1 .CB are 6.5, 12.0 and 7.05 angstroms respectively, and the 
distances of "Zipper 2 " -F.ALA1 .CB-B.ALA4. CB-F.ALA3. CB-B.ALA6. CB are 6.6, 
11.79, 6.96 angstroms respectively. 5b shows the far vdW surface. 5c shows 
the violet colored ABEF Chains of Figure 3. 5d shows HBs: A/E.ALA5. 0- 
B/F.GLY2.N, A/E.ALA5.N-B/F. GLY2. 0, A/E.ALA3. 0-B/F.ALA4.N, A/E.ALA3.N- 
B/F.ALA4.0, A/E.ALA1.0-B/F.ALA6.N. 
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-752.075, -285.005, 3576.69) and its Hessian matrix H = 



/ 66.4497 -10.415 -27.1845 \ 

-10.415 82.2093 54.6557 

-27.1845 54.6557 203.929 

380.664 -4.02618 -159.578 

-4.02618 369.586 -25.5081 

\ -159.578 -25.5081 1048.02 / 



which is a positive definite matrix with eigenvalues (1085.02, 372.093, 341.157, 
230.049, 61.2695, 61.2695). Then Algorithm 1 hybridized with simulated an- 
nealing global optimal search (in order to bring local optimal solutions to jump 
out of local traps, replacing the discrete gradient local optimal search method in 
Algorithm 1 of [fl~4l by the Algorithm 1 of this paper) is executed and the optimal 
solution (3.027, 4.954, 3.856; 1.679, 1.777, 5.011) for Eq. 3.5 is got. 

Similarly, for Eq. 3.7, we take the initial solution x =(4.714, 4.878, 8.881; 
4.170, 1.360, 14.292) and calculate its gradient c/ =(-46.8782, 133.142, 397.798, 
-401.192, -182.604, 3436.46) and its Hessian matrix H = 



/ 66.1163 -6.0961 -18.2092 \ 

-6.0961 81.3119 51.7918 

-18.2092 51.7918 218.677 

339.302 -2.9188 -85.8315 

-2.9188 361.487 -2.99786 

\ -85.8315 -2.99786 1034.38 / 



which is a positive definite matrix with eigenvalues (1044.83, 361.8, 328.538, 
238.159, 63.973, 63.973). The optimal solution (7.412, 4.760, 4.624; 5.887, 
1.451, 5.757) for Eq. 3.7 is got. 

We set (3.027, 4.954, 3.856; 1.679, 1.777, 5.01 1) as the coordinates of E.ALA6.CB 
and E.ALA4.CB, (7.412, 4.760, 4.624; 5.887, 1.451, 5.757) as the coordinates of 
F.ALA1.CB and F.ALA3.CB, and taking the average value we get 

E(F) = A(B) + (1.0335, 1.0823, 0.9723). (3.10) 

By Eq. 3.10 we can get very close vdW contacts between A (B) chains and E(F) 
chains (Figure 5(b)), and other chains of Model 1 can be got by Eqs. 3.2~3.3 and 
Eq. 3.11: 

I (J) = A{E) - (9.666, 0, 0). (3.11) 
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Figure 6: Left: The close vdW surface contacts between Chains A(B) and E(F) of 
Model 1 after solving Eqs. 3.5 and 3.7; Right: The constructed Model 1 with 10 
chains. The dashed lines denote HBs. 

The initial structure of Model 1 illuminated in Figure 6(a) ~6(b) is not the op- 
timal structure with the lowest total potential energy. The initial structure also has 
no hydrogen atoms (so no hydrogen bonds existed) and water molecules added. 
For each Chain, the C-terminal and N-terminal atoms also have problems. Clearly 
there are a lot of close/bad contacts between /3-strand atoms as illuminated in 
Figure 6(a)~6(b). We used the ff03 force field of AMBER 11, in a neutral pH 
environment. The amyloid fibrils were surrounded with a 8 angstroms layer of 
TIP3PBOX water molecules using the XLEaP module of AMBER 11. 1944 wa- 
ters and 408 hydrogen atoms were added for Model 1 byt he XLEaP module, 
solvated amyloid fibril was inimized by the method. 

The LJ potential energy of atoms' vdW interactions is just a part of the to- 
tal potential energy of a protein, and by observations from Model 1 computed 
by Eqs. 3.10, 3.2~3.3 and 3.11 we can see that the contacts between /3-strand 
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atoms and /3-sheet atoms are too close/bad. Thus, we need to relax Model 1 com- 
puted. The relaxation is done in the use of local search LBFGS Quasi-Newton 
method (lbfgs_memory_depth=3) within AMBER 11 ll23l . The relaxed/optimized 
Model 1 is illuminated in Figure 6(c). Seeing Figure 6(d) compared with Figure 
6(b), we may know the vdW interactions between the two /3-sheets are very re- 
laxed/optimized now. Figure 6(c) shows the Model 1 of optimal molecular struc- 
ture for prion AGAAAAGA amyloid fibrils. 

4 Conclusion 

In a (macro)molecular system, a robust energy minimization is very necessarily 
required. Mathematical optimization minimization methods find a place to apply 
in these systems. Because in physics the (macro)molecular system usually is not 
a simple two-body problem of system, local search optimization methods are very 
useful in the applications to the (macro)molecular system. On anther sense, when 
a protein is unstable, noncrystalline or insoluble and very difficult to detect its 3D 
structure by the expensive and costly NMR and X-ray, theoretical mathematical or 
physical computational method can be used to produce the 3D structure of the pro- 
tein. Moreover, even the X-ray crystallography finds the X-ray final structure of a 
protein, we still need refinements using theoretical protocols in order to produce 
a better structure. The theoretical computational method - an improved LBFGS 
Quasi-Newtonian mathematical optimization method - presented in this paper and 
other mathematical optimization methods mentioned in this paper should be very 
useful in the protein molecular modeling research field. 

This paper also shows the effectiveness of the improved LBFGS mathemat- 
ical optimization method presented. Prion AGAAAAGA amyloid fibrils have 
not much structural information. This paper presents some bioinformatics on the 
molecular structures of prion AGAAAAGA amyloid fibrils in the sense of theo- 
retical emphasis. The structures may be helpful in the advance in the biochemical 
knowledge of prion protein misfolding or instability and in the future applications 
for therapeutic agent design. 



Acknowledgments 

This research was supported by the National Natural Science Foundation of China 
(No. 1 1 171 180), the National Basic Research Program of China (2010CB732501), 
and a Victorian Life Sciences Computation Initiative (VLSCI) grant numbered 



16 



VR0063 & 488 on its Peak Computing Facility at the University of Melbourne, 
an initiative of the Victorian Government (Australia). The first author appreciates 
the financial support by Professors Changyu Wang and Yiju Wang for his visits 
to Qufu Normal University (China), and by Professor Xiangsun Zhang for his 
visits to Chinese Academy of Sciences (Beijing); the first author also appreciates 
Professor Adil M. Bagirov (University of Ballarat) for his instructions to imple- 
ment the Wolfe-type search. This paper is dedicated to Professor Xiangsun Zhang 
(Academia Sinica, Beijing) on the occasion of his 70th birthday. 

References 

[1] Chiang PK, Lam MA, Luo Y (2008) The many faces of amyloid beta in 
Alzheimer's disease. Curr Mol Med 8(6): 580-584. 

[2] Irvine GB, El-Agnaf OM, Shankar GM, Walsh DM (2008) Protein aggre- 
gation in the brain: the molecular basis for Alzheimer's and Parkinson's 
diseases. Mol med Cambridge Mass 14(7-8): 451-64. 

[3] Ferreira ST, Vieira MN, De Felice FG (2007) Soluble protein oligomers as 
emerging toxins in Alzheimer's and other amyloid diseases. IUBMB life 
59(4-5): 332-45. 

[4] Nature Editorial (2001) More than just mad cow disease. Nature Struct Biol 
8(4): 281. doi: 10.1038/86132 

[5] Truant R, Atwal RS, Desmond C, Munsie L, Tran T (2008) Huntington's 
disease: revisiting the aggregation hypothesis in polyglutamine neurode- 
generative diseases. The FEBS J 275(17): 4252-4262. 

[6] Weydt P, La Spada AR (2006) Targeting protein aggregation in 
neurodegeneration-lessons from polyglutamine disorders. Expert opinion 
on therapeutic targets 10(4): 505-513. 

[7] Nelson R, Sawaya MR, Balbimie M, Madsen A, Riekel C, Grothe R, Eisen- 
berg D (2005) Structure of the cross-beta spine of amyloid-like fibrils. Na- 
ture 435(7043): 773-778. 

[8] Sawaya MR, Sambashivan S, Nelson R, Ivanova MI, Sievers SA, Apos- 
tol MI, Thompson MJ, Balbirnie M, Wiltzius JJ, McFarlane HT, Madsen 
A, Riekel C, Eisenberg D (2007) Atomic structures of amyloid cross-beta 
spines reveal varied steric zippers. Nature 447(7143): 453-457. 



17 



[9] Zhang JP (2011) The Lennard- Jones potential minimization problem for 
prion AGAAAAGA amyloid fibril molecular modeling. arXiv: 1 1 06. 15841 

[10] Kolossvry I, Bowers KJ (2010) Global optimization of additive potential 
energy functions: predicting binary Lennard-Jones clusters. Phys Rev E 
Stat Nonlin Soft Matter Phys 82(5 Pt 2): 05671 1. 

[11] Sicher M, Mohr S, Goedecker S (2010) Efficient moves for global geometry 
optimization methods and their application to binary systems. J Chem Phys 
134(4): 44-106. 

[12] Strodel B, Lee JW, Whittleston CS, Wales DJ (2010) Transmembrane 
structures for Alzheimer's A/3(l-42) oligomers. J Am Chem Soc 132(38): 
13300-13312. 

[13] Ye T, Xu R, Huang W (201 1) Global optimization of binary Lennard-Jones 
clusters using three perturbation operators. J Chem Inf Model 51(3): 572- 
577. 

[14] Zhang JP, Sun J, Wu CZ (2011) Optimal atomic-resolution structures of 
prion AGAAAAGA amyloid fibrils. J Theor Biol 279(1): 17-28. 

[15] Zhang JP, Gao DY, Yearwood J. (2011) A novel canonical dual computa- 
tional approach for prion AGAAAAGA amyloid fibril molecular modeling. 
J Theor Biol 284(1): 149-157. 

[16] Grosso A, Locatelli M, Schoen F (2009) Solving molecular distance geom- 
etry problems by global optimization algorithms. Comput Optim Appl 43: 
23-37. 

[17] More JJ, Wu ZJ (1997) Global continuation for distance geometry prob- 
lems. SIAM J Optim 7: 814-836. 

[18] Zou ZH, Bird RH, Schnabel RB (1997) A stochastic/perturbation global 
optimization algorithm for distance geometry problems. J Glob Optim 11: 
91-105. 

[19] Brown DR (2000) Prion protein peptides: optimal toxicity and peptide 
blockade of toxicity. Mol Cell Neurosci 15: 66-78. 

[20] van der Spoel D, Lindahl E, Hess B, van Buuren AR, Apol E, Meulenhoff 
PJ, Tieleman DP, Sijbers ALTM, Feenstra KA, van Drunen R, Berendsen 
HJC (2010) Gromacs User Manual version 4.5.4, www.gromacs.org. 



18 



[21] Snyman JA (2005) Practical Mathematical Optimization: An Introduction 
to Basic Optimization Theory and Classical and New Gradient-Based Al- 
gorithms. Springer Publishing, ISBN 0-387-24348-8. 

[22] Chu, J, Trout BL, Brooks BR (2003) A super-linear minimization scheme 
for the nudged elastic band method. J Chem Phys 119: 12708-12717. 

[23] Case DA, Darden TA, Cheatham TE, Simmerling III CL, Wang J, Duke 
RE, Luo R, Walker RC, Zhang W, Merz KM, Roberts BP, Wang B, Hayik 
S, Roitberg A, Seabra G, Kolossvary I, Wong KF, Paesani F, Vanicek J, Liu 
J, Wu X, Brozell SR, Steinbrecher T, Gohlke H, Cai Q, Ye X, Wang J, Hsieh 
M-J, Cui G, Roe DR, Mathews DH, Seetin MG, Sagui C, Babin V, Luchko 
T, Gusarov S, Kovalenko A, Kollman PA (2010) AMBER 11, University of 
California, San Francisco. 

[24] Liu DC, Nocedal J (1989) On the limited memory method for large scale 
optimization. Math Programming B 45: 503-528. 

[25] Nocedal J, Morales J (2000) Automatic preconditioning by limited memory 
quasi-Newton updating. SIAM J Opt 10: 1079-1096. 

[26] Byrd RH, Nocedal J, Zhu C (1995) Towards a discrete Newton method with 
memory for large-scale optimization, Optimization Technology Center Re- 
port OTC-95-1, EEC Department, North-western University. 

[27] Byrd RH, Lu P, Nocedal J (1995) A limited memory algorithm for bound 
constrained optimization. SIAM J Scientif Statistic Comput 16: 1 190— 
1208. 

[28] Zhu C, Byrd RH, Nocedal J (1997) L-BFGS-B: Algorithm 778: L-BFGS- 
B, FORTRAN routines for large scale bound constrained optimization. 
ACM Trans Math Softw 23: 550-560. 

[29] Dennis JE, Robert JR, Schnabel B (1996) Numerical Methods for Uncon- 
strained Optimization and Nonlinear Equations, SIAM. 

[30] Nocedal J, Wright SJ (1999) Numerical Optimization, Springer- Verlag, 
ISBN 0-387-98793-2 & Numerical Optimization (2nd ed.), Berlin, New 
York: Springer- Verlag, ISBN 978-0-387-30303-1 

[31] Berndt E, Hall B, Hall R, Hausman J (1974) Estimation and Inference in 
Nonlinear Structural Models. Annals of Economic and Social Measurement 
3: 653-665. 



19 



[32] Luenberger DG, Ye YY (2008) Linear and nonlinear programming (3rd 
edition), International Series in Operations Research & Management Sci- 
ence 116, New York: Springer, ISBN 978-0-387-74502-2. 

[33] Wang YJ, Xiu NH (2012) Nonlinear Optimization Theory and Methods, 
Science Press, Beijing, ISBN 9787030330284. 

[34] Powell MJD (1984) Nonconvex minimization calculations and the conju- 
gate gradient method. In: Griffiths D.F. ed., Numerical Analysis, Lecture 
Notes in Mathematics 1066, Springer Verlag, Berlin: 122-141. 

[35] Mascarenhas WF (2004) The BFGS method with exact line searches fails 
for non-convex objective functions. Math Program Ser A 99: 49-61. 

[36] Dai YH (2002) Convergence properties of the BFGS Algorithm. SIAM J 
Optiml3(3): 693-701. 

[37] Li DH, Fukushima M (2001) A modified BFGS method and its global con- 
vergence in nonconvex minimization. J Comput Math Appl 129: 15-35. 

[38] Xiao YH, Wei ZX, Wang ZG (2008) A limited memory BFGS-type method 
for large-scale unconstrained optimization. J Comput Math Appl 56(4): 
1001-1009. 

[39] Yang YT, Xu CX (2007) A compact limited memory method for large scale 
unconstrained optimization. Eur J Oper Res 180(1): 48-56. 

[40] Hou YT, Wang YJ (2012) The modified limited memory BFGS method 
for large-scale optimization, MSc degree thesis, Qufu Normal University, 
China. 

[41] Grippo L, Lamparello F, Lucidi S (1986) A nonmonotone line search tech- 
nique for Newton's method. SIAM J Numer Anal 23: 707-716. 

[42] Han JY, Liu GH (1997) Global convergence analysis of a new nonmonotone 
BFGS algorithm on convex objective Functions. Compu Optim Appl 7: 
277-289. 

[43] Bagirov AM, Karaszen B, Sezer M (2008) Discrete gradient method: 
derivative-free method for nonsmooth optimization. J Optim Theor Appl 
137(2): 317-334. 

[44] Dolan ED, Mor JJ (2002) Benchmarking optimization software with per- 
formance profiles. Math Program Ser A 91: 201-213. 



20 



[45] Byrd RH, Nocedal J (1989) A tool for the analysis of quasi-Newton meth- 
ods with application to unconstrained minimization. SIAM J Numer Anal 
26: 727-739. 

[46] Byrd RH, Nocedal J, Yuan YX (1987) Global convergence of a class of 
quasi-Newton methods on convex problems. SIAM J Numer Anal 24: 
1171-1189. 

[47] Dennis JE, More JJ (1977) Quasi-Newton method, motivation and theory. 
SIAM Rev 19: 46-89. 

[48] Nocedal J (1980) Updating quasi-Newton matrices with limited storage. 
Math Comput 35: 773-782. 

[49] Gilbert JC, Lemarichal C (1989) Some numerical experiments with vari- 
able storage quasi-Newton algorithms. Math Programming 45: 407-435. 

[50] Al-Baali M (1999) Improved Hessian approximations for the limited mem- 
ory BFGS method. Numer. Algorithms 22: 99-112. 

[51] Byrd RH, Nocedal J, Schnabel RB (1994) Representations of quasi-Newton 
matrices and their use in limited memory methods. Math Programming 63: 
129-156. 

[52] Yuan GL, Wei ZX, Wu YL (2010) Modified limited memory BFGS method 
with nonmonotone line search for unconstrained optimization. J Korean 
Math Soc 47: 767-788. 

[53] Alper T, Cramp W, Haig D, Clarke M (1967) Does the agent of scrapie 
replicate without nucleic acid?. Nature 214(5090): 764-766. 

[54] Griffith J (1967) Self-replication and scrapie. Nature 215(5105): 1043- 
1044. 

[55] Brown DR (2001) Microglia and prion disease. Microsc Res Tech 54: 71- 
80. 

[56] Brown DR, Herms J, Kretzschmar HA (1994) Mouse cortical cells lacking 
cellular PrP survive in culture with a neurotoxic PrP fragment. Neuroreport 
5: 2057-2060. 

[57] Cappai R, Collins SJ (2004) Structural biology of prions, in Prions A Chal- 
lenge for Science, Medicine and the Public Health System (Rabenau H.F., 
Cinatl J., Doerr H.W. (eds)), Basel, Karger 11: 14-32. 



21 



[58] Chabry J, Caughey B, Chesebro B (1998) Specific inhibition of in vitro 
formation of protease-resistant prion protein by synthetic peptides. J Biol 
Chem 273(21): 13203-13207. 

[59] Cheng HM, Tsai TWT, Huang WYC, Lee HK, Lian HY, Chou FC, Mou 
Y Chu J, Chan JC (2011) Steric zipper formed by hydrophobic peptide 
fragment of Syrian hamster prion protein. Biochem 50(32): 6815-6823. 

[60] Gasset M, Baldwin MA, Lloyd DH, Gabriel JM, Holtzman DM, Cohen 
F, Fletterick R, Prusiner SB (1992) Predicted alpha- helical regions of the 
prion protein when synthesized as peptides form amyloid. Proc Natl Acad 
Sci USA 89(22): 10940-10944. 

[61] Haigh CL, Edwards K, Brown DR (2005) Copper binding is the governing 
determinant of prion protein turnover. Mol Cell Neurosci 30(2): 186-196. 

[62] Harrison CF, Lawson VA, Coleman BM, Kim YS, Masters CL, Cappai 
R, Barnham KJ, Hill AF (2010) Conservation of a glycine-rich region in 
the prion protein is required for uptake of prion infectivity. J Biol Chem 
285(26): 20213-20223. 

[63] Holscher C, Delius H, Burkle A (1998) Overexpression of nonconvertible 
PrP c deltal 14-121 in scrapie-infected mouse neuroblastoma cells leads 
to trans-dominant inhibition of wild-type PrP Sc accumulation. J Virol 72: 
1153-1159. 

[64] Jobling MF, Huang X, Stewart LR, Barnham KJ, Curtain C, Volitakis I, Pe- 
rugini M, White AR, Cherny RA, Masters CL, Barrow CJ, Collins SJ, Bush 
AI, Cappai R (2001) Copper and Zinc binding modulates the aggregation 
and neurotoxic properties of the prion peptide PrP 106-126. Biochem 40: 
8073-8084. 

[65] Jobling MF, Stewart LR, White AR, McLean C, Friedhuber A, Maher F, 
Beyreuther K, Masters CL, Barrow CJ, Collins SJ, Cappai R (1999) The 
hydrophobic core sequence modulates the neurotoxic and secondary struc- 
ture properties of the prion peptide 106-126. J Neurochem 73: 1557-1565. 

[66] Jones EM, Wu B, Surewicz K, Nadaud PS, Helmus JJ, Chen S, Jaroniec CP, 
Surewicz WK (2011) Structural polymorphism in amyloids: new insights 
from studies with Y145Stop prion protein fibrils. J Biol Chem 286(49): 
42777-42784. 



22 



[67] Kourie JI (2001) Mechanisms of prion-induced modifications in membrane 
transport properties: implications for signal transduction and neurotoxicity. 
Chem Biol Interact 138(1): 1-26. 

[68] Kourie JI, Kenna BL, Tew D, Jobling MF, Curtain CC, Masters CL, Barn- 
ham KJ, Cappai R (2003) Copper modulation of ion channels of PrP[106- 
126] mutant prion peptide fragments. J Membr Biol 193(1): 35-45. 

[69] Kuwata K, Matumoto T, Cheng H, Nagayama K, James TL, Roder H 
(2003) NMR-detected hydrogen exchange and molecular dynamics simu- 
lations provide structural insight into fibril formation of prion protein frag- 
ment 106-126. Proc Natl Acad Sci USA 100(25): 14790-14795. 

[70] Laganowsky A, Liu C, Sawaya MR, Whitelegge JP, Park J, Zhao M, Pen- 
salfini A, Soriaga AB, Landau M, Teng PK, Cascio D, Glabe C, Eisenberg 
D (2012) Atomic view of a toxic amyloid small oligomer. Science 335: 
1228-1231. 

[71] Lee SW, Mou Y, Lin SY, Chou FC, Tseng WH, Chen CH, Lu CY, Yu SS, 
Chan JC (2008) Steric zipper of the amyloid fibrils formed by residues 109- 
122 of the Syrian hamster prion protein. J Mol Biol 378(5): 1 142-1 154. 

[72] Ma BY, Nussinov R (2002) Molecular dynamics simulations of alanine rich 
/3-sheet oligomers: insight into amyloid formation. Protein Sci 11: 2335- 
2350. 

[73] Norstrom EM, Mastrianni JA (2005) The AGAAAAGA palindrome in PrP 
is required to generate a productive PrP Sc -PrP c complex that leads to prion 
propagation. J Biol Chem 280: 27236-27243. 

[74] Sasaki K, Gaikwad J, Hashiguchi S, Kubota T, Sugimura K, Kremer W, 
Kalbitzer HR, Akasaka K (2008) Reversible monomer-oligomer transition 
in human prion protein. Prion 2(3): 1 18-122. 

[75] Wagoner VA (2010) Computer simulation studies of self-assembly of fib- 
ril forming peptides with an Intermediate resolution protein model. PhD 
thesis, North Carolina State University, Raleigh, North Carolina. 

[76] Wagoner VA, Cheon M, Chang I, Hall CK (2011) Computer simulation 
study of amyloid fibril formation by palindromic sequences in prion pep- 
tides. Proteins Struct Funct Bioinf 79(7): 2132-2145. 

[77] Wegner C, Romer A, Schmalzbauer R, Lorenz H, Windl O, Kretzschmar 
HA (2002) Mutant prion protein acquires resistance to protease in mouse 
neuroblastoma cells. J Gen Virol 83: 1237-1245. 



23 



[78] Zanuy D, Ma B, Nussinov R (2003) Short peptide amyloid organization: 
stabilities and conformations of the islet amyloid peptide NFGAIL. Bio- 
physJ 84(3): 1884-1894. 



[79] Zhang JP (2011) Optimal molecular structures of prion AGAAAAGA 
palindrome amyloid fibrils formatted by simulated annealing. J Mol Model 
17(1): 173-179. 

[80] Zhang ZQ, Chen H, Lai LH (2007) Identification of amyloid fibril-forming 
segments based on structure and residue-based statistical potential. Bioinf 
23: 2218-2225. 

Supplementary Materials [40|| 
Mathematical Assumptions (A1)-(A4): 

(Al) the objective function / : R™ — > R is twice continuously differentiable, 
i.e. feC 2 ; 

(A2) / has Lipschitz continuous gradients and Hessians, i.e., there exist con- 
stants L, M 2 > such that 

\\g{x) - g{y)\\ < L\\x - y\\, Vi,?/Gl n , 

\\G(x) - G(x*)\\ < M 2 ||x-x*||, 

where G denotes the Hessian and x* is a stationary point (i.e. x fc converges to x* 
where g(x*) = 0); 

(A3) G(x*) is positive definite; and 

(A4) the level set of / 

Q = {x e R n | f(x) < f(x )}, 
is bounded, where x E R™ is an initial point. 

Mathematical Proof of the Algorithm Convergence and Conver- 
gent Rate: 

In order to establish the convergence for Algorithm 1, we represent Eq. 2.3 as 



b l n k b l ilk b l 
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where B k = H k 1 ,s l = x l+1 - x h y\ = y t + 7,s,, and B k k ™ +1 = B for all k. 



Theorem 1 Let Assumptions (A1)^(A4) hold, and {x k } denotes the sequence 
generated by Algorithm , 

stantt G [0, 1) such that 



generated by Algorithm 1. Then lim inf = 0; moreover, there exits a con 



f(x k )-f(x*)<t k [f(x -f(x*)i 

which is to say that x k converges to the minimum R-linearly. 

Proof. Denote by M 1 an upper bound of \\gk\\ on f2 and by y* k = y k + 7fcSfc, we 
can prove that 

\\yff < (L±M 1 l 
Vk s k 7 

where 7 is a sufficiently small positive number. By this formula and taking the 
trace operation in both sides of Eq. 5. 1, we get 

tr(B k+1 ) = tr(B^)- E f£ + E ^ 

l =k -m+l 1 l=k-m+l Ul 

< tr(B )-if£ + ilf£ (5.2) 

< tr(5„)-E^ + «(Hl) 

1=0 1 

< M 3 , 

where M 3 is a positive constant. We also take the determinant and by Eq. 5.2, we 
have 



k-m+l\ TT Vi ' 

11 sfB 

■k-fh+1 1 

k «t ......... (5.3) 



det(B k+1 ) = det{Bt m+1 ) Y[ ^ 

l=k-m+l 1 

det{B,)Y[^>det{B,){^ k ^ 
1=0 1 



> M 4 , 

where M 4 is a positive constant. Combining Eqs. 5.2 and 5.3, we obtain a con- 
stant 5 > such that cos(6 k ) > 5, where cos(6> fe ) - s * BkSk 



\\sk\\\\B k s k \ 



Because g(x) is Lipschitz continuous and by the Wolfe-type line search con- 
dition Eq. 2.2, we have 

(o- 2 - l)g k dk < (9k+i - 9k) T d k < A fc L||4|| 2 , 
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which implies 



Using the Wolfe-type line search condition Eq. 2.1 and Eq. 5.4, we have 

fk+i < fk + <Ti*k9%d k <f k + a 1 (CT2 - 1} ( £f 



\\d k \\ 

ill II, „ 

Then using the monotonicity and the bound of {f k } on fi, we obtain 



/fe + <7l^||<y 2 COS 2 fc . 



||^fc|| 2 COS 2 fc < 00. 

fc=l 

Combining cos 6^ > 5, we get 

lim inf \\gk\\ = 0. 



Under the Assumptions (A2)~(A3), there is a neighbourhood £/(x*) of x* 
such that for all x e £/ (x*), 

1 1 (oj) | j > ||<?(x) — <?(x*)|| > m||x — x*\\, (5.5) 

and for all d e iT, 

d T G(x)d > m||d|| 2 , (5.6) 

where m is a positive constant got according to the Assumptions. Hence, Eqs. 
5.5~5.6 hold with x = x^ for all k sufficiently large. 

Since 

yf s k = (Vk + lkS k ) T s k > yls k 

> -(1 - a 2 )gls k 

= (1 - <7 2 )||#fc||||Sfc|| cos 9 k , 

and y* k T s k < (L + Mi) 2 ||s fe || 2 , if denoting cci = ( L ^^ 2 , we can obtain 

||sfc|| > ,j , *l\i \\9k\\ cos6 k = ai\\g k \\ cos6 k . (5.7) 

Combining Eq. 5.7 and the Wolfe-type line search condition Eq. 2.1, we get 

fk+i-f* < fk~f* + ^i9k s k 

= /k-/*-<7l||0k||||Sfc||cOS0* 

< fk-f*-(Jiai\\g k \\ 2 cos 2 e k . 
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Using Eq. 5.5, then Eq. 5.7 implies for all k sufficiently large, 

fk+i-f* < fk - f* - o-iailbfcll 2 cos 2 e k , 5g . 

By Taylor's expansion, then there exists a positive number M 5 > such that 

f k -r<M 5 \\x k -x*\\ 2 

for all k sufficiently large. 

Let a 2 = £1 §f ! , indeed Eq. 5.8 implies that 

fk+i-f* < fk- f* ~ (Jiaim 2 ||a; fc - x*|| 2 cos 2 ^ 
< fk- r-a 2 (f k - f*) cos 2 9 k 
= (l-a 2 cos 2 e k )(f k -f*) 

for all sufficiently large k. 

Hence, there is a constant t = 1 — a 2 cos 2 9 k e [0, 1) such that 

f(x k )-f(x*)<t k [f(x )-f(x*)\. (5.9) 



By Eq. 5.6, we have 

1 



fk~ f* > 2 m W X k - x *\\ 2 ' 



which, together with Eq. 5.9, implies 

\\x k — x*\\ <t 



fc/2r 2(/(x ) -f(x*) 1/2 
m 



This is to say the sequence {x k } is R- linearly convergent. The proof of this theo- 
rem is complete. ■ 
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